home
***
CD-ROM
|
disk
|
FTP
|
other
***
search
/
The Datafile PD-CD 1 Issue 2
/
PDCD-1 - Issue 02.iso
/
_utilities
/
utilities
/
001
/
meschach
/
!Meschach
/
c
/
arnoldi
next >
Wrap
Text File
|
1994-01-13
|
5KB
|
188 lines
/**************************************************************************
**
** Copyright (C) 1993 David E. Steward & Zbigniew Leyk, all rights reserved.
**
** Meschach Library
**
** This Meschach Library is provided "as is" without any express
** or implied warranty of any kind with respect to this software.
** In particular the authors shall not be liable for any direct,
** indirect, special, incidental or consequential damages arising
** in any way from use of the software.
**
** Everyone is granted permission to copy, modify and redistribute this
** Meschach Library, provided:
** 1. All copies contain this copyright notice.
** 2. All modified copies shall carry a notice stating who
** made the last modification and the date of such modification.
** 3. No charge is made for this software or works derived from it.
** This clause shall not be construed as constraining other software
** distributed on the same medium as this software, nor is a
** distribution fee considered a charge.
**
***************************************************************************/
/*
Arnoldi method for finding eigenvalues of large non-symmetric
matrices
*/
#include <stdio.h>
#include <math.h>
#include "matrix.h"
#include "matrix2.h"
#include "sparse.h"
static char rcsid[] = "$Id: arnoldi.c,v 1.3 1994/01/13 05:45:40 des Exp $";
/* arnoldi -- an implementation of the Arnoldi method */
MAT *arnoldi(A,A_param,x0,m,h_rem,Q,H)
VEC *(*A)();
void *A_param;
VEC *x0;
int m;
Real *h_rem;
MAT *Q, *H;
{
static VEC *v=VNULL, *u=VNULL, *r=VNULL, *s=VNULL, *tmp=VNULL;
int i;
Real h_val;
if ( ! A || ! Q || ! x0 )
error(E_NULL,"arnoldi");
if ( m <= 0 )
error(E_BOUNDS,"arnoldi");
if ( Q->n != x0->dim || Q->m != m )
error(E_SIZES,"arnoldi");
m_zero(Q);
H = m_resize(H,m,m);
m_zero(H);
u = v_resize(u,x0->dim);
v = v_resize(v,x0->dim);
r = v_resize(r,m);
s = v_resize(s,m);
tmp = v_resize(tmp,x0->dim);
MEM_STAT_REG(u,TYPE_VEC);
MEM_STAT_REG(v,TYPE_VEC);
MEM_STAT_REG(r,TYPE_VEC);
MEM_STAT_REG(s,TYPE_VEC);
MEM_STAT_REG(tmp,TYPE_VEC);
sv_mlt(1.0/v_norm2(x0),x0,v);
for ( i = 0; i < m; i++ )
{
set_row(Q,i,v);
u = (*A)(A_param,v,u);
r = mv_mlt(Q,u,r);
tmp = vm_mlt(Q,r,tmp);
v_sub(u,tmp,u);
h_val = v_norm2(u);
/* if u == 0 then we have an exact subspace */
if ( h_val == 0.0 )
{
*h_rem = h_val;
return H;
}
/* iterative refinement -- ensures near orthogonality */
do {
s = mv_mlt(Q,u,s);
tmp = vm_mlt(Q,s,tmp);
v_sub(u,tmp,u);
v_add(r,s,r);
} while ( v_norm2(s) > 0.1*(h_val = v_norm2(u)) );
/* now that u is nearly orthogonal to Q, update H */
set_col(H,i,r);
if ( i == m-1 )
{
*h_rem = h_val;
continue;
}
/* H->me[i+1][i] = h_val; */
m_set_val(H,i+1,i,h_val);
sv_mlt(1.0/h_val,u,v);
}
return H;
}
/* sp_arnoldi -- uses arnoldi() with an explicit representation of A */
MAT *sp_arnoldi(A,x0,m,h_rem,Q,H)
SPMAT *A;
VEC *x0;
int m;
Real *h_rem;
MAT *Q, *H;
{ return arnoldi(sp_mv_mlt,A,x0,m,h_rem,Q,H); }
/* gmres -- generalised minimum residual algorithm of Saad & Schultz
SIAM J. Sci. Stat. Comp. v.7, pp.856--869 (1986)
-- y is overwritten with the solution */
VEC *gmres(A,A_param,m,Q,R,b,tol,x)
VEC *(*A)();
void *A_param;
VEC *b, *x;
int m;
MAT *Q, *R;
double tol;
{
static VEC *v=VNULL, *u=VNULL, *r=VNULL, *tmp=VNULL, *rhs=VNULL;
static VEC *diag=VNULL, *beta=VNULL;
int i;
Real h_val, norm_b;
if ( ! A || ! Q || ! b || ! R )
error(E_NULL,"gmres");
if ( m <= 0 )
error(E_BOUNDS,"gmres");
if ( Q->n != b->dim || Q->m != m )
error(E_SIZES,"gmres");
x = v_copy(b,x);
m_zero(Q);
R = m_resize(R,m+1,m);
m_zero(R);
u = v_resize(u,x->dim);
v = v_resize(v,x->dim);
tmp = v_resize(tmp,x->dim);
rhs = v_resize(rhs,m+1);
MEM_STAT_REG(u,TYPE_VEC);
MEM_STAT_REG(v,TYPE_VEC);
MEM_STAT_REG(r,TYPE_VEC);
MEM_STAT_REG(tmp,TYPE_VEC);
MEM_STAT_REG(rhs,TYPE_VEC);
norm_b = v_norm2(x);
if ( norm_b == 0.0 )
error(E_RANGE,"gmres");
sv_mlt(1.0/norm_b,x,v);
for ( i = 0; i < m; i++ )
{
set_row(Q,i,v);
tracecatch(u = (*A)(A_param,v,u),"gmres");
r = mv_mlt(Q,u,r);
tmp = vm_mlt(Q,r,tmp);
v_sub(u,tmp,u);
h_val = v_norm2(u);
set_col(R,i,r);
R->me[i+1][i] = h_val;
sv_mlt(1.0/h_val,u,v);
}
/* use i x i submatrix of R */
R = m_resize(R,i+1,i);
rhs = v_resize(rhs,i+1);
v_zero(rhs);
rhs->ve[0] = norm_b;
tmp = v_resize(tmp,i);
diag = v_resize(diag,i+1);
beta = v_resize(beta,i+1);
MEM_STAT_REG(beta,TYPE_VEC);
MEM_STAT_REG(diag,TYPE_VEC);
QRfactor(R,diag /* ,beta */);
tmp = QRsolve(R,diag, /* beta, */ rhs,tmp);
v_resize(tmp,m);
vm_mlt(Q,tmp,x);
return x;
}